Structure of cold nuclear matter at subnuclear densities by Quantum Molecular 

Dynamics 
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Structure of cold nuclear matter at subnuclear densities for the proton fraction x = 0.5, 0.3 and 
0.1 is investigated by quantum molecular dynamics (QMD) simulations. We demonstrate that the 
phases with slablike and rodlike nuclei, etc. can be formed dynamically from hot uniform nuclear 
matter without any assumptions on nuclear shape, and also systematically analyze the structure of 
cold matter using two-point correlation functions and Minkowski functionals. In our simulations, 
we also observe intermediate phases, which has complicated nuclear shapes. It is found out that 
these phases can be characterized as those with negative Euler characteristic. Our result implies the 
existence of these kinds of phases in addition to the simple "pasta" phases in neutron star crusts 
and supernova inner cores. 

In addition, we investigate the properties of the effective QMD interaction used in the present 
work to examine the validity of our results. The resultant energy per nucleon e„ of the pure neutron 
matter, the proton chemical /ip"' in pure neutron matter and the nuclear surface tension -Egm-f are 
generally reasonable in comparison with other nuclear interactions. 

PACS numbers: 21.65.+f,26.50.+x,26.60.+c,61.20. Ja 



I. INTRODUCTION 

For the past several decades since the discovery of pul- 
sars, many authors have investigated the properties of 
dense matter which exist inside neutron stars and su- 
pernova cores (see, e.g., Refs. These objects 
have been shown that they consist of a variety of mate- 
rial phases whose physical properties reflect in many as- 
trophysical phenomena of these objects. Especially, the 
properties of nuclear matter under extreme conditions, 
which is one of the essential topics for understanding the 
mechanism of collapse-driven supernovae 0|, the struc- 
ture of neutron star crusts |^ and its relating phenomena, 
have been studied actively. This subject is also interest- 
ing as one of the fundamental problems of the complex 
fluids of nucleons. 

At subnuclear densities, nuclear matter exhibits the 
coexistence of a liquid phase with a gas phase due to the 
internucleon interaction which has an attractive part. At 
sufficiently low temperatures relevant to neutron star in- 
teriors, and sufficiently below the normal nuclear density, 
long-range Coulomb interactions make the system divide 
periodically into gas and spherical liquid drops, adding a 
crystalline property to the liquid-gas coexistence. 

In the density region where nuclei are about to melt 
into uniform matter, it is expected that the energetically 
favorable configuration of the mixed phase possesses in- 
teresting spatial structures such as rodlike and slablike 
nuclei and rodlike and spherical bubbles, etc., which are 
referred to as nuclear "pasta" . This picture was origi- 
nally proposed by Ravenhall et al. 6] and Hashimoto 
et al. independently. Their predictions both based 
on free energy calculations with liquid drop models as- 



suming some specific nuclear shapes. These works clarify 
that the most energetically stable nuclear shape is deter- 
mined by a subtle balance between the nuclear surface 
and Coulomb energies. Detailed aspects of equilibrium 
phase diagrams, such as a series of nuclear shapes which 
can be realized as the energetically most favorable and 
the density range corresponding to the phases with non- 
spherical nuclei, vary with nuclear models 0. However, 
the realization of the "pasta" phases as energy minimum 
states can be seen in a wide range of nuclear models and 
the phase diagrams possess an universal basic feature 
that, with increasing density, the shape of the nuclear 
matter region changes like sphere —^ cylinder — > slab — > 
cylindrical hole spherical hole uniform ^| . This 
feature is also reproduced by the Thomas-Fermi calcula- 
tions by several groups [ll|, [13, UM ■ 

The phases with these exotic nuclear structures, if they 
were realized in neutron star crusts or supernova cores, 
bring about many astrophysical consequences. As for 
those in neutron star phenomena, it is interesting to note 
the relevance of nonspherical nuclei to pulsar glitches and 
cooling of neutron stars. Although the question whether 
the mechanism of pulsar glitches is depicted by vortex 
pinning model or star quake model has yet to be set- 
tled completely, the existence of nonspherical nuclei in 
neutron star matter (NSM) have significant effects in 
both cases. As for the former, while the force needed 
to pin vortices has yet to be clarified completely even 
for a bcc lattice of spherical nuclei mainly due to the 
uncertain properties of impurities and defects |14| . the 
effect of spatial structure of normal nuclear matter on 
vortex dynamics cannot be ignored. As for the latter, 
the existence of "pasta" phases with slablike and rodlike 
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nuclei would change the elastic properties of inner crust 
matter from those of crystalline solid to those of liquid 
crystal as indicated by Pethick and Potekhin , which 
results in significant decrease of the maximum elastic en- 
ergy that can be stored in the inner crust. The presence 
of nonspherical nuclei would also accelerate the cooling 
of the corresponding region of neutron stars by opening 
semileptonic weak processes which are unlikely to occur 
for spherical nuclei 

"Pasta" phases in supernova matter (SNM) are ex- 
pected to affect the neutrino transport and hydrodynam- 
ics in supernova cores. Let us first note that the neutrino 
wavelengths, typically of order 20 fm, are comparable 
to or even greater than the internuclear spacing, lead- 
ing to diffractive effects on the neutrino elastic scatter- 
ing off such a periodic spatial structure of nuclear matter 
|g. These effects, induced by the internuclear Coulombic 
correlations, would reduce the scattering rates and hence 
the lepton fraction Yl- For the bcc lattice of spheri- 
cal nuclei, such a reduction was examined by Horowitz 
[T^ by calculating the associated static structure factor. 
It is also noteworthy that nonspherical nuclei and bub- 
bles are elongated in specific direction. In such direction, 
the neutrino scattering processes are no longer coherent, 
in contrast to the case of roughly spherical nuclei whose 
finiteness in any direction yields constructive interference 
in the scattering. The final point to be mentioned is that 
the changes in the nuclear shape are accompanied by dis- 
continuities in the adiabatic index, denoting how hard the 
equation of state of the material is. These discontinuities 
may influence the cor e hy drodynamics during the initial 
phase of the collapse [Tjl ■ 

Though the properties of "pasta" phases in equilibrium 
state have been investigated actively, the formation and 
the melting processes of these phases have not been dis- 
cussed except for some limited cases which are based on 
perturbative approaches 0,^3- It is important to adopt 
a microscopic and dynamical approach which allows ar- 
bitrary nuclear structures in order to understand these 
processes of nonspherical nuclei. At finite temperatures, 
it is considered that not only nuclear surface becomes 
obscure but also nuclei of various shapes may coexist. 
Therefore, it is necessary to incorporate density fluctua- 
tions without any assumptions on nuclear shape to inves- 
tigate the properties of "pasta" phases at finite tempera- 
tures. Although the works done by Williams and Koonin 
[Tl| and Lassaut et al. ^3 do not assume nuclear struc- 
ture, they can not incorporate fluctuations of nuclcon 
distributions in a satisfying level because they are based 
on the Thomas-Fermi calculation, which is one-body ap- 
proximation. In addition, only a single structure is con- 
tained in the simulation box in these works, there are 
thus possibilities that nuclear shape is strongly affected 
by boundary effect and some structures are prohibited 
implicitly. 

In the present work, we study the structure of cold 
dense matter at subnuclear densities in the framework 
of quantum molecular dynamics (QMD) |l8j |. which is 



one of the molecular dynamics (MD) approaches for nu- 
cleon many-body systems (see, e.g., Ref. jfJl for re- 
view). MD for nucleons including QMD, which is a mi- 
croscopic and dynamical method without any assump- 
tions on nuclear structure, is suitable for incorporating 
fluctuations of particle distributions. Previously, we have 
reported the first results of our study on nuclear "pasta" 
by QMD, which demonstrated that the "pasta" phases 
can be formed in a dynamical wa y fo r matter with pro- 
ton fractions x = 0.3 and 0.5 |23|. In this full pa- 
per, we present new results for astrophysically interesting 
neutron-rich matter of x = 0.1 in addition to the cases 
of a; = 0.3 and 0.5 reported before. 

The plan of this paper is as follows. In Section II, we 
describe the framework of the QMD model used in the 
present study. We then show the results of our simula- 
tions in Section III and analyze the structure of matter 
obtained by the simulations using two-point correlation 
functions and Minkowski functionals in Section IV. In 
Section V, we investigate the properties of the effective 
nuclear interaction used in this work in order to examine 
the validity of our results in terms of nuclear forces. As- 
trophysical discussions are given in Section VI. Summary 
and conclusions are presented in Section VII. 



II. QUANTUM MOLECULAR DYNAMICS 

We have various types of molecular dynamics meth- 
ods for nucleons including representative ones such as 
fermionic molecular dynamics (FMD) j23|, antisym- 
metrized molecular dynamics (AMD) [231 and QMD, etc. 
In the present work, we choose QMD among them bal- 
ancing between calculation cost and accuracy. The typ- 
ical length scale / of inter-structure is Z ~ 10 fm and 
the density region of interest is just below the normal 
nuclear density po — 0.165 fm~'^. The required nucleon 
number N in order to reproduce n unit structures in the 
simulation box is about N ^ po{nl)'^ (for slabs). It is 
thus desirable that we prepare nucleons of order 10000 
if we try to reduce boundary effects down to a satisfac- 
tory level by reproducing several several unit structures 
in the box. While it is a hard task to treat such a large 
system with, for example, FMD and AMD whose calcula- 
tion costs scale as ~ N'^, it is feasible to do it with QMD 
whose calculation costs scale as ^ N"^. This difference 
comes from summations in the Slater determinants in 
the trial wave functions of the former models. In QMD, 
on the other hand, the total A^- nucleon wave function |<I>) 
is assumed to be a direct product of single-nucleon states 

1$) = ® |(/)2) • • • ® . (1) 

The single-nucleon state is represented by a Gaussian 
wave packet: 
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where Ri(t) and Pi(i) are the centers of position and 
momentum of the packet j, respectively, and L is a pa- 
rameter related to the extension of the wave packet in 
the coordinate space. 

It is also noted that we mainly focus on the macro- 
scopic structures; the exchange effect would not be so 
important for them. This can be seen by comparing the 
typical values of the exchange energy for the macroscopic 
scale and of the energy difference between two successive 
phases with nonspherical nuclei. Suppose there are two 
identical nucleons, i ^ 1 and 2, bound in different nuclei 
each other. The exchange energy between these particles 
is calculated as an exchange integral: 

K = J U{ri - r2) (pi(ri)(^i(r2)¥32(r2)<^2(i"i) dridr2 , 

(3) 

where U is the potential energy. An asymptotic form of 
the wave function is given by 



exp {-Kir) 



(4) 



with r = |ri — and Ki = j-\/2mEi , {i = 1,2), where 
Ei is the binding energy and m is the nucleon mass. The 
exchange integral reads 



K ~ exp + K2)R] 5 X lO^^McV 



(5) 



for the internuclear distance i? ~ 10 fm and Ei ~ 8 MeV, 
which is extremely smaller than the typical energy differ- 
ence per nucleon between the "pasta" phases of order 0.1 
keV (for NSM, see Fig. 4 in Ref. 0) - 10 keV (for SNM, 



see Fig. 4 in Ref. ^3)- Therefore, it is expected that 
QMD, which is less elaborate in treating the exchange 
effect, is not bad approximation for investigating the nu- 
clear "pasta". Consequently, QMD has the advantages 
over the other models in the present study. In the fu- 
ture, we will have to confirm the validity of the results 
obtained by QMD using other more elaborate model such 
as AMD or FMD to treat the exchange effect more pre- 
cicely. However, this problem is beyond the scope of the 
present work. 



A. Model Hamiltonian 

To simulate nuclear matter at subnuclear densities 
within the framework of QMD, we use a QMD model 
Hamiltonian developed by Maruyama et al. which 
is constructed so as to reproduce bulk properties of nu- 
clear matter and properties of finite nuclei. This model 
Hamiltonian consists of the following six terms 



H — T + Vpauli + Vskyimc + Vsy 



- Vmd + Vco 



ilomb 



, (6) 



where T is the kinetic energy, Vpauii is the Pauli potential 
introduced to reproduce the Pauli principle effectively, 
Vskyrmc is thc Skyrme potential which consists of an at- 
tractive two-body term and a repulsive three-body term, 
Vsym is the symmetry potential, Vmu is the momentum- 
dependent potential introduced as two Fock terms of the 
Yukawa interaction and Vcouiomb is the Coulomb poten- 
tial. The expressions of these terms are given as 
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(9) 
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(11) 
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where pij means the overlap between the single-nucleon as 

densities, pi{r) and Pj{r), for z-th and j-th nucleons given . 

= d r p,{r) pj{r) , (13) 
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CTi is the nucleon spin and is the isospin 
{Ti = 1/2 for protons and —1/2 for neutrons) and 

Cp, 90, Po, a, 13, T, Cs, Cci\ Ccx , Mi, and L are 
model parameters determined to reproduce the proper- 
ties of the ground states of the finite nuclei, especially 
heavier ones, and the saturation properties of nuclear 
matter '2^ . A parameter set used in this work is shown 
in Table m The single- nucleon densities Pi{r) and Pi{r) 
are given by 



with 



\Mr)\' 



(27rL)3/2 



1 



exp 



(r-R.)' 



(27ri)3/2 

(r-R.) 



2L 



exp 



2L 



(14) 
(15) 

(16) 



The modified width L in pi(r) is introduced in the three- 
body term of Skyrme interaction [Eq. @] to incorporate 
the effect of the repulsive density-dependent term. 



TABLE L Effective interaction parameter set 
(incompressibility K=280 MeV; medium EOS model in Ref. 

Hi) 



Cp (MeV) 


207 


po (MeV/c) 


120 


go (fm) 


1.644 


a (MeV) 


-92.86 


(3 (MeV) 


169.28 


T 


1.33333 


Cs (MeV) 


25.0 


Cii' (MeV) 


-258.54 


dx' (MeV) 


375.6 


Ml (fm-i) 


2.35 


M2 (fm-i) 


0.4 


L (fm^) 


2.1 



We adopt QMD equations of motion with friction 
terms to simulate the dynamical relaxation: 



(17) 



where the friction coefhcients and are positive def- 
inite, which determine the relaxation time scale. The 
relaxation scheme given by Eqs. (|17|l is referred to as the 
steepest descent method and it leads to the continuous 
decrease in H as 



~dt 



Ri 



P. 



f dn 



.dn 



< 



(18) 



Even though it is recognized that this method is not ef- 
ficient, it is expected that the dynamics given by Eqs. 
((T7jl with S,R,S,p <^ 1 deviates slightly in a short period 
from the physically grounded dynamics given by QMD 
equations of motion without the friction terms [equations 
without the second terms in the right-hand sides of Eqs. 
(|17|l ]. which we would like to respect. 



III. QMD SIMULATIONS OF COLD MATTER 
AT SUBNUCLEAR DENSITIES 

A. QMD Simulations for x = 0.5 and 0.3 

We have performed QMD simulations of an infinite 
{n,p,e) system with fixed proton fractions a; = 0.5 and 
0.3 for various nucleon densities p [the density region is 
(0.05 - 1.0) Po]- We set 2048 nucleons (1372 nucleons 
in some cases) contained in a cubic box which is im- 
posed by the periodic boundary condition. Throughout 
this paper, the numbers of the protons (neutrons) with 
up-spin and with down-spin are equal. The relativistic 
degenerate electrons which ensure the charge neutrality 
are regarded as a uniform background and the Coulomb 
interaction is calculated by the Ewald method taking ac- 
count of the Gaussian charge distribution of each wave- 
packet (see Appendix^. This method enables us to effi- 
ciently sum up contributions of long-range interactions in 
a system with periodic boundary conditions. For nuclear 
interaction, we use the effective Hamiltonian developed 
by Maruyama et al. (medium EOS model) m whose 
expressions are given in the last Section. 

We first prepare an uniform hot nucleon gas at A:bT ~ 
20 MeV as an initial condition equilibrated for ~ 500 — 
2000 fm/c in advance. In order to realize the ground 
state of matter, we then cool it down slowly for 0(10'^ — 
10**) fm/c, keeping the nucleon density constant with the 
frictional relaxation method [Eqs. H17|l ]. etc. |3| until 
the temperature gets ~ 0.1 MeV or less. Note that no 
artificial fluctuations are given in the simulation. 

The QMD equations of motion with the friction terms 
given by Eqs. ((T7jl are solved using the fourth-order Gear 
predictor-corrector method in conjunction with multi- 
ple time step algorithm |26j. Integration time steps At 
are set to be adaptive in the range of At < 0.1 — 0.2 
fm/c depending on the degree of convergence. At each 
step, the correcting operation is iterated until the er- 
ror of position Ar and the relative error of momen- 
tum Ap/p become smaller than 10~®, where Ar and 
Ap/p are estimated as the maximum values of correc- 
tion among all particles. We mainly use PCs(Pentium 
III) equipped with MDGRAPE-2, which accelerates cal- 
culations of momentum-independent forces including the 
long-range Coulomb force. 

Shown in Figs. ^ and |21 are the resultant nucleon dis- 
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FIG. 1: (Color) The nucleon distributions of typical phases with simple structures of cold matter at x — 0.5; (a) sphere 
phase, O.lpo (ibox = 43.65 fm, N = 1372); (b) cylinder phase, 0.225po (ibox = 38.07 fm, N = 2048); (c) slab phase, 0.4po 
(Lbox = 31-42 fm, TV = 2048); (d) cylindrical hole phase, 0.5po (Lbox = 29.17 fm, iV = 2048) and (e) spherical hole phase, 0.6po 
(ibox = 27.45 fm, TV — 2048), where Lbox is the box size. The red particles represent protons and the green ones represent 
neutrons. 




FIG. 2; (Color) Same as Fig. [T] at x = 0.3; (a) sphere phase, O.lpo (i^box = 49.88 fm, TV = 2048); (b) cylinder phase, 
0.18po (I'box = 41.01 fm, TV = 2048); (c) slab phase, 0.35po (ibox = 32.85 fm, TV = 2048); (d) cylindrical hole phase, 0.5po 
(ibox = 29.17 fm, TV = 2048) and (e) spherical hole phase, 0.55/9o (^box = 28.26 fm, TV = 2048). The red particles represent 
protons and the green ones represent neutrons. 
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(b) x=0.3 
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^ <Q BP: sphere CH: cylindrical hole 
C : cylinder SH: spherical hole 
S : slab ( , ): intermediate phase 

FIG. 3: Phase diagrams of cold matter at a; = 0.5 (a) and x — 
0.3 (b). Matter is unstable against phase separation in the 
density region shown as kt < 0, where kt is the isothermal 
compressibility. The symbols SP, C, S, CH and SH stand 
for nuclear shapes, i.e., sphere, cylinder, slab, cylindrical hole 
and spherical hole, respectively. The parentheses (A,B) show 
intermediate phases between A-phase and B-phase suggested 
in this work. They have complicated structures different from 
those of both A-phase and B-phase. Simulations have been 
carried out at densities denoted by small circles. 



shown in Figs.l^a) and (b) for x = 0.5 and 0.3, respec- 
tively. As can be seen from these figures, the obtained 
phase diagrams basically reproduce the sequence of the 
energetically favored nuclear shapes predicted by simple 
discussions (7j vifhich only take account of the Coulomb 
and surface effects; this prediction is that the nuclear 
shape changes like sphere — s- cylinder — s- slab cylin- 
drical hole — > spherical hole — > uniform, with increasing 
density. Comparing Figs. 3 (a) and (b), we can see that 
the phase diagram shifts towards the lower density side 
with decreasing x, which is due to the tendency that the 
saturation density is lowered as the neutron excess in- 
creases. It is remarkable that the density dependence of 
the nuclear shape, except for spherical nuclei and bub- 
bles, is quite sensitive and phases with intermediate nu- 
clear shapes which are not simple as shown in Figs. 1 and 
2 are observed in two density regions: one is between the 
cylinder phase and the slab phase, the other is between 
the slab phase and the cylindrical hole phase. We note 
that these phases are different from coexistence phases 
with nuclei of simple shapes, which will be referred to as 
"intermediate phases" . 



tributions of cold matter at a; = 0.5 and 0.3, respectively. 
We can see from these figures that the phases with rod- 
like and slablike nuclei, cylindrical and spherical bubbles, 
in addition to the phase with spherical nuclei are repro- 
duced in both the cases of a; = 0.5 and 0.3. We here would 
like to mention the reasons of discrepancies between the 
present result and the result obtained by Maruyama et 
al. which says "the nuclear shape may not have these 
simple symmetries" . One of the most crucial reasons 
seems to be the difference in treatment of the Coulomb 
interaction. In the present simulation, we calculate the 
long range Coulomb interaction in a consistent way us- 
ing the Ewald method. For the system of interest where 
the Thomas-Fermi screening length is comparable to or 
larger than the size of nuclei, this treatment is more ad- 
equate than that which inctroduces an artificial cutoff 
distance as in Ref. The other crucial reason is the 

difference in the relaxation time scales r fm/c; we set 
T ^ 0(10^ — 10"*) in the present work, but Maruyama 
et al. set t ~ severalxlO^ fm/c j23|. In our simulation, 
we can reproduce the bubble-phases [see (d) and (e) of 
Figs. Hand 12] with r ~ 10"^ fm/c and the nucleus-phases 
[see (b) and (c) of Figs, d and |21 with r ~ 0(10^^) fm/c. 
However, the matter in the density region correspond- 
ing to a nucleus-phase is quenched in an amorphous-like 
state when r < 10^ fm/c. In the present work, we take 
T much larger than typical time scale rth ~ O(IOO) fm/c 
for nucleons to thermally diffuse in the distance of Z ^ 10 
fm at /9 ~ po and fceT" — 1 MeV. This temperature is 
lower than the typical value of the liquid-gas phase tran- 
sition temperature in the density region of interest, it is 
thus considered that our results are thermally relaxed in 
a satisfying level. 

Phase diagrams of matter in the ground state are 



B. QMD Simulations for x = 0.1 

We have also performed QMD simulations of matter 
with proton fraction a; = 0.1 as a more realistic condition 
for the neutron star matter. In this case, we have to deal 
with a larger system than in the cases of x = 0.5 and 
0.3 because enough number of protons for reproducing 
several nuclei in a simulation box are required to obtain 
significant results; protons play an important role in gen- 
erating the long-range order due to their electric charge. 
We have investigated the neutron-rich matter at a; = 0.1 
with 10976 nucleons, in which 1098 protons and 9878 
neutrons are contained. Following basically the same 
procedure that used for the cases of x = 0.5 and 0.3 (see 
Section Fill Al for detail), we tried to obtain the ground- 
state matter. However, in the present case, we quickly 
relax from the initial state at fceT ^ 20 MeV to the 
state at ksT ^ 10 MeV, at which matter is still uniform, 
with a Nose-Hoover-like thermostat |2^ |2^ which will 
be discussed in another paper [2^. After the relaxation 
at fceT - 10 MeV for - 4000 - 7000 fm/c, we then cool 
down the system with a relaxation time scale r ~ 0(10*^) 
fm/c using the QMD equations of motion with friction 
terms (|17|l . These simulations are performed by Fujitsu 
VPP 5000 equipped with MDGRAPE-2. 

Some resultant nucleon distributions are shown in Figs. 
0] and O which correspond to the sphere phase and the 
cylinder phase, respectively. As can be seen in Fig. 0] 
dripped neutrons spread over the whole region in the 
simulation box, which lead to smaller density contrast 
compared with that for the cases of a; = 0.5 and 0.3 de- 
picted in Figs. Hand [3 respectively. 

The results obtained for x = 0.1 are summarized in the 
phase diagram shown in Fig. |^ A striking feature is that 
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FIG. 4: (Color) The nucleon distribution of sphere phase in 
cold matter at x = 0.1. The nucleon density p and the size 
ibox of the simulation box are p — 0.075po and Lbox ~ 96.08 
fm. The red particles represent protons and the green ones 
represent neutrons. 



the intermediate phase (C,U)]. However, in the present 
neutron-rich case, "pasta" phase with slablike nuclei can- 
not be obtained as far as we have investigated, which 
will be discussed at the end of the next section. It is 
also noted that the density at which matter turns into 
uniform is lower than those in the cases of x = 0.5 and 
0.3, which is consistent with the tendency that matter 
becomes the more neutron-rich, the saturation density 
becomes the lower. 

x=0.1 



SP 


c 




(C,U) 


uniform 


0.1 0.2 


0.3 0.4 0.5 


0.6 0.7 0.8 0.9 1.0 



SP: sphere 
C : cylinder 

U : uniform ( , ): intermediate ptiase 

FIG. 6: Phase diagram of cold matter at x = 0.1. Mat- 
ter is unstable against phase separation in the density region 
shown as kt < 0, where kt is the isothermal compressibility. 
The symbols SP, C and U stand for shapes of nuclear matter 
region, i.e., sphere, cylinder and uniform, respectively. The 
density at which the area-averaged mean curvature of nuclear 
surface is zero is denoted by (H) — 0. However, slab phase is 
not observed in our results even at such a density. The paren- 
theses (A,B) show intermediate phases between A-phase and 
B-phase. They have complicated structures different from 
those of both A-phase and B-phase. Simulations have been 
carried out at densities denoted by small circles. 



IV. ANALYSIS OF THE STRUCTURE OF 
MATTER 

A. Two-point correlation functions 

To analyze the spatial distribution of nucleons, we cal- 
culate two-point correlation function for nucleon den- 
sity field p'*-* (i = N,p,n; where N stands for nucleons). 
^ii is here defined as 



FIG. 5; (Color) The proton distribution of cylinder phase 
in cold matter at x = 0.1. The nucleon density p and the size 
ibox of the simulation box are p = 0.2po and Lbox ~ 69.29 fm. 
Neutrons which spread over the whole space are not depicted 
in this figure. 



the wide density region from ~ 0.25po to ^ 0.525po is oc- 
cupied by an intermediate phase. The structure of matter 
seems to change rather continuously from that consists 
of branching rodlike nuclei connected to each other [ob- 
tained in the lower density region of the intermediate 
phase denoted by (C,U)] to that consists of branching 
bubbles connected to each other [higher density region of 



1- fdn,^ I d^x 5,(x)5,;(x + r) (19) 



(S,{x)St (x + r))x,o. 



(20) 



where (• • ■)x,nr denotes an average over the position x 
and the direction of r, and (5i(x) is the fluctuation of the 
density field p^'^ (x) given by 



,_ pW(x)-p(0 

= == 



with 



,(■0 



V 



(21) 



(22) 
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FIG. 7: Two-point correlation functions of density fluctua- 
tions calculated for x = 0.5. The solid lines show the two- 
point correlation function for nucleon density distributions; 
the dotted lines, that for proton density distributions; the 
dashed lines, that for neutron density distributions. 



We construct the nucleon density distribution p^*^(x) 
from a data set of the centers of position of the nucle- 
ons by the following procedure. We first set 64'^ (for 
X — 0.5 and 0.3) or 128^ (for x = 0.1) grid points in 
the simulation box and then distribute particle num- 
bers on each grid point using the cloud-in-ccU method 
(see, e.g., Ref. H^)- Next, we carry out the smoothing 
procedure in the discrete Fourier space with a Gaussian 
smoothing function corresponding to the distribution of 
the wave packet given by Eq. The density distribu- 
tions p^^n^) — SfcLi '^fc''l'^fc(^)p7 where n[,*^ = or 1 
projects on particle type i, in the discrete real space can 
be obtained by the inverse Fourier transformation. The 
Fourier transformations are performed using the FFT al- 
gorithm. 

The resultant two-point correlation functions £,NNir), 
fpp(r) and Cnnif) at various densities below at which 
matter becomes uniform at zero temperature for x — 0.5, 
0.3 and 0.1 are plotted in Figs. [7| |H1 and 1^1 respectively. 
We can see the general tendency, which is common for 
the different values of x, that the amplitude of £,ii{r) de- 
creases with increasing the density. It is noted that even 
though the change in the amplitude of (r) is quite no- 
ticeable, the smallest zero-point r = ro of i,u{r) takes 
similar values at various densities especially for x — 0.5 
and 0.3. This feature means that the typical length scales 
of the nuclear structures, i.e., the internuclear distance 
and the nuclear radius, keep comparable at subnuclear 
densities from ~ O.lpo to pm, which is consistent with 
the results obtained by the previous works (see, e.g., Refs. 
H H m m). This behavior just below pn, will be dis- 
cussed further concerning a problem about the properties 
of the transition to uniform matter. 

We also note that a strong attractive force acting be- 
tween a proton and a neutron leads to the good agree- 
ment of the zero-points of and ^„„ even for x — 0.3 
and 0.1 as well as for x = 0.5 although the zero-point 
''o of £,nn is ~ 0.3 fm (^ 0.5 fm) larger than that of 






FIG. 8: Two-point correlation functions of density fluctua- 
tions calculated for x = 0.3; (a) for nucleon density distri- 
butions; (b) for proton density distributions; (c) for neutron 
density distributions. 



for X = 0.3 (0.1) at each density. This shows that the 
phase of the density fluctuations of protons and neutrons 
correlate so strongly with each other at zero temperature 
that they almost coincide. 

As can be seen by comparing S^NNir) for different val- 
ues of the proton fraction x, the amplitude decreases and 
the value of the smallest zero-point vq increases with de- 
creasing X. This behavior means that, as matter becomes 
more neutron-rich, not only the nucleon density distribu- 
tion gets smoother but also the spatial structure becomes 
larger. 

Let us then examine for each value of the proton 
fraction. For symmetric matter (x = 0.5), protons and 
neutrons are equivalent except for the mass difference 
and the Coulomb interaction. Therefore, the two-point 
correlation functions £,nni Cpp a-i^d Crm are almost the 
same at larger values of r > tq. At smaller values of 
r, ^pp is slightly smaller than i^„„ because the repulsive 
Coulomb interaction among protons tends to reduce the 
proton density inhomogeneity especially in the smaller 
scale. For asymmetric matter {x = 0.3 and 0.1), on the 
contrary, the amplitude of ^„„ is much smaller than that 
of ^pp due to the dripped neutrons which distribute rather 
uniformly outside the nuclei. 



B. Transition to uniform matter 

Let us here examine the properties of the transition 
from the phase with spherical bubbles to uniform matter 
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FIG. 9; The same as Fig. |Hlfor x = 0.1. 



for X — 0.5. For this purpose, two-point correlation func- 
tion of the nucleon density fluctuation is useful. In Fig. 
1101 we thus plot the two-point correlation function of the 
nucleon density fluctuation £,NN{f) for several densities 
around the melting density p^. To compute S,nn{t)i we 
use a 1372-nucleon system cooled down until the temper- 
ature gets ^0.05 MeV by QMD equations of motion H17() 
with friction terms. 

For miiform phase, (,nn{t) should be zero except for 
the contribution of short-range correlation. The behavior 
of ^ArAr(r) shows that Pin hes between 0.7po and 0.725po 
at a; = 0.5 [see also Fig.|2Ia)] above which long-range cor- 
relation disappears. It is noted that the smallest value 
of r — tq at which ^NN{r) = keeps around 8 fm even 
at densities just below p^. This means that, for x = 0.5, 
the phase with spherical bubbles whose radii are around 
ro suddenly disappear rather than they shrink gradually 
and the system turns into uniform with increasing the 
density because the quantity tq nearly corresponds to 
the half wave-length of the inhomogeneous density pro- 
file. The discontinuous change in the density profile indi- 
cates that the transition between the phase with spheri- 
cal bubbles and the uniform phase is of first order. This 
conclusion is also obtained in the previous calculations 
for which the spatial structure of the nuclear matter re- 
gion and/ or the shape of the density profile are assumed 
0, 0, 0, 0, 1^. In the present work, we have con- 
firmed the first order nature by QMD simulations with- 
out these assumptions as several authors have done so 
without the assumptions by the Thomas-Fermi approxi- 
mation [Hill. 

For the cases of a; = 0.3 and 0.1, we could not see the 
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FIG. 10: Two-point correlation function ^NN{r) of the nu- 
cleon density fluctuation. ^]vjv(r) is calculated for x = 0.5 
and for densities around pm from 0.625po to 0.75po. 

significant sign of the first order nature of the transition 
between the mixed phase and the uniform phase because 
the amplitude of £,NNii") is quite small just below p^. 
Further study is necessary to determine the properties of 
the transition for theses cases of asymmetric matter. 

C. Minkowski functionals 

To extract the morphological characteristics of the nu- 
clear shape changes and the intermediate phases, we in- 
troduce the Minkowski functionals (see, e.g., Ref. [s^ 
and references therein; a concise review is provided by 
Ref. '33^) as geometrical and topological measures of the 
nuclear surface. Let us consider a homogeneous body 
K ^ TZ m the c?-dimensional Eucledian space, where TZ 
is the class of such bodies. Morphological measures are 
defined as functionals (p : TZ ^ R which satisfy the fol- 
lowing three general properties: 

(1) Motion invariance : The functional is indepen- 
dent of the position and the direction of the body, 
i.e., 

<p(X) = ^{gK) , (23) 
where g denotes any translations and rotations. 

(2) Additivity : The functional of the union of two 
bodies should behave like a volume. The contribu- 
tion of the overlapping region should be subtracted 
, i.e., 

(p(Xi U K2) = ifiKi) + <p(if2) - <p(ifi n K2) , (24) 
where Ki, K2 E TZ- 

(3) Continuity : If the body is approximated with pix- 
els, the functional of the approximate body con- 
verges to that of the original body when the pixels 
get smaller , i.e., 

lim Lp{Kn) = (p{K) as lim ^ K , (25) 

n — >-00 71 — >oo 
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where K is a convex body and {Kn} is a sequence 
of convex bodies. 

Hadwiger's theorem in integral geometry states that 
there are just d + 1 independent functionals which sat- 
isfy the above properties; they are known as Minkowski 
functionals. In three dimensional space, four Minkowski 
functionals are related to the volume, the surface area, 
the integral mean curvature and the Euler characteristic. 

In classifying nuclear shapes including those of the in- 
termediate phases obtained in our simulations, the in- 
tegral mean curvature and the Euler characteristic are 
useful, which will be discussed later. Both are described 
by surface integrals of the following local quantities: the 
mean curvature H = {ki + K2)/2 and the Gaussian cur- 
vature G = K1K2, i.e., Jgj^HdA and x = 5^ lax^'^^^ 
where ki and K2 are the principal curvatures and dA is 
the area element of the surface of the body K. The Eu- 
ler characteristic x is a purely topological quantity and 
is given by 

X = (number of isolated regions) — (number of tunnels) 
-I- (number of cavities). (26) 

Thus X > for the sphere and the spherical hole phases 
and the coexistence phase of spheres and cylinders, and 
X = for the other ideal "pasta" phases, i.e. the cylin- 
der, the slab and the cylindrical hole phases which con- 
sist of infinitely long rods, infinitely extending slabs and 
infinitely long cylindrical holes, respectively. We intro- 
duce the area-averaged mean curvature, {H) = ^/ HdA, 
and the Euler characteristic density, x/^j normalized 
quantities, where V is the volume of the whole space. 



1. Minkowski functionals for x = 0.5 and 0.3 

We calculate the normalized Minkowski functionals, 
i.e., the volume fraction u, the surface area density A/V, 
the area-averaged mean curvature (H) and the Euler 
characteristic density x/^ for ^ = 0-5 and 0.3 by the fol- 
lowing procedure. As described in Subsection IIV Al we 
first construct proton and nucleon density distributions 



x-0.5 



x-0.3 



Ap) 



|0..(r)pand p{^) = ELiIM^ 



where n'l!'' = or 1 is the isospin projection on the 



Ap) _ 

proton state. We set a threshold proton density Pp^th 
and then calculate /(pp,th) = V'(pp,th)M(Pp,th), where 
l^(Pp,th) and A{pp,th) are the volume and the surface area 
of the regions in which p'^P^(r) > Pp,th- We find out the 
value pp,th = P*,th where 



/(Pp.th) = and define 



the regions in which ^'■^''(r) > p* as nuclear regions. 
For spherical nuclei, for example, p* corresponds to 
a point of inflection of a radial density distribution. In 
the most phase-separating region, the values of p^j^ dis- 
tribute in the range of about 0.07 — 0.09 fm~'^ in the both 
cases of X = 0.5 and 0.3, where p^^^ is the threshold nu- 
cleon density corresponds to p* t^- We then calculate u, 
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FIG. 11: Density dependence of the volume fraction u and 
the surface area density A/V of cold matter at a; = 0.5 (a and 
c) and X — 0.3 (b and d). The crosses show the results for 
pth = Pth and the open triangles and squares show the results 
for pth = Pth ~ 0.05po and p'^ -I- 0.05po, respectively. 
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FIG. 12: Density dependence of the area-averaged mean cur- 
vature (H) and the Euler characteristic density x/^ of cold 
matter at a; = 0.5 (a and c) and x = 0.3 (b and d). 



A, J HdA and x for the identified nuclear surface. We 
evaluate u by counting the number of pixels at which 
p(P^(r) is higher than p* ^i^, A by the triangle decomposi- 
tion method, / HdA by the algorithm shown in Ref. [s^l 
in conjunction with a calibration by correction of surface 
area, and x by the algorithm of Ref. and by that 
of counting deficit angles 0|, which are confirmed that 
both of them give the same results. 
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We have plotted the resultant p dependence of u, A/V , 
{H) and x/V" for the isodcnsity surface of pth ~ Pth 
in Figs. ^3 and ^1 In addition to the values of u, 
A/V and (H) for the isodcnsity surface of pth = Pth' 
have also investigated those for the isodcnsity surfaces of 
Pth — Pth^O.OSpo to examine the extent of the uncertain- 
tics of these quantities which stem from the arbitrariness 
in the definition of the nuclear surface. As shown in Fig. 
111! these uncertainties are at most ~ 0.1 for u and ~ 0.25 
fm~^ for A/V. For (H), we could not observe remark- 
able differences from the values for pth = Pth (they were 
smaller than 0.015 fm~^). We could not see these kinds of 
uncertainties in x/^i except for the densities near below 

Pni- 

As shown in Fig. ^2 the volume fraction u of the nu- 
clear regions increases almost monotonically under pm in 
both cases of a; = 0.5 and 0.3. This feature reflects the 
incompressible nature of nuclear matter. It is interest- 
ing to see the density dependence of the nuclear surface 
density A/V because this quantity is directly related to 
the surface energy density, which is one of the key fac- 
tors in determining the nuclear shape. Figs. Ilir cl and 
lllf dl show that, as the nucleon density increases, A/V 
increases at a nearly constant rate until p ~ O.Spoj a-nd 
then its increasing rate becomes rather smaller around 
the density region of the slab phase, and finally it be- 
gins to decrease in the density region of the cylindrical 
hole phase or the spherical hole phase. This general be- 
havior can be understood from the density dependence 
of the surface energy density obtained by simple argu- 
ments, which only allow for the nuclear surface and the 
Coulomb effects (see, e.g., Refs. 0,l23)- 

We also plot u and A/V for the nucleon density dis- 
tribution p(r) as functions of the threshold density pth 
evaluated at various values of p. The results for x — 0.5 
and 0.3 are shown in Figs. 1131 and 1141 respectively. Fea- 
tures of nuclear shape changes can be seen in the behav- 
ior of the curves of A/V. Peaks in the higher pth region 
are attributed to nucleons in the nuclear matter regions 
and broad bumps in the lower pth region around 0.05 - 
O.lpo observed for x = 0.3 are due to the dripped neu- 
trons outside nuclei. A/V in the intermediate pth region 
in which its slope is nearly constant mainly comes from 
contribution of nuclear surfaces. As the nucleon density 
p increases, the higher pth peak becomes more clear and 
the position of the center of the peak finally coincides 
with p in the uniform phase. This feature shows that 
the nuclear matter regions become more uniform with 
increasing the density. When the dispersion of the inter- 
nucleon distance is small, large surface area caused by 
the Gaussian density distribution of each nucleon can be 
picked up with a single value of pth- As can be seen in 
Fig. 1141 the lower pth bump, in turn, disappears with in- 
creasing p. This is because, as p increases, the dripped 
neutron gas becomes more inhomogeneous and tends to 
distribute close to the nuclear surface leading to a lower 
proton fraction in the nuclear matter regions. It is also 
noted that, as the nucleon density increases, the slope in 
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FIG. 13: Volume fraction u (upper panel) and surface area 
density A/V (lower panel) as functions of the threshold den- 
sity pth calculated for x = 0.5 and various nucleon densities 
P- 



the intermediate pth region changes from negative to pos- 
itive at the density corresponding to the phase of slablike 
nuclei (0.4po for x = 0.5 and 0.35po for x = 0.3), which 
is consistent with what is expected from the sign of (H) 
for the nuclear surface [see Figs. I12r a) and ll2f b)]. 

Let us then focus on {H) and x/F to classify the nu- 
clear shape. The behavior of (H) shows that it decreases 
almost monotonically from positive to negative with in- 
creasing p until the matter turns into uniform. The den- 
sities corresponding to (H) ~ are about 0.4 and 0.35po 
for X — 0.5 and 0.3, respectively; these values are consis- 
tent with the density regions of the phase with slablike 
nuclei (see Fig. As mentioned previously, x/^ is ac- 
tually positive in the density regions corresponding to 
the phases with spherical nuclei, coexistence of spherical 
and cylindrical nuclei, and spherical holes because of the 
existence of isolated regions. As for those corresponding 
to the phases with cylindrical nuclei, planar nuclei and 
cylindrical holes, x/^ — 0- The fact that the values of 
x/y are not exactly zero for nucleon distributions shown 
as the slab phase in Figs. and [3 reflects the imperfec- 
tion of these "slabs" , which is due to the small nuclear 
parts connecting the neighboring slabs. However, we can 
say that the behaviors of x/^ plotted in Figs. I12r c) and 
I12f d) show that x/V is negative in the density region of 
the intermediate phases, even if we take into account the 
imperfection of the obtained nuclear shapes and the un- 
certainties of the definition of the nuclear surface. This 
means that the intermediate phases consist of nuclear 
surfaces which are saddle-like at each point on average 
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FIG. 14: The same as Fig. [SI for x = 0.3. 

and they consist of each highly connected nuclear and 
gas regions due to a lot of tunnels [see Eq. H26|) ]. Using 
the quantities (H) and x/^i the sequence of the nuclear 
shapes with increasing the density can be described as 
follows: ((H) > 0, x/V > 0) ^ {{H) > 0, x/V = 0) 
^ {{H) > 0, x/V < 0) ^ ((F> = 0, x/V = 0) 
^ ((H) < 0, x/V < 0) ^ ((H) < 0, x/V = 0) ^ 
((H) < 0, x/V > 0) ^ uniform. 

Let us now consider the discrepancy from the results 
of previous works which do not assume nuclear structure 
[ill lia | ; the intermediate phases can not be seen in these 
works. We can give following two reasons for the discrep- 
ancy. 

(1) These previous calculations are based on the 
Thomas-Fermi approximation which can not suf- 
ficiently incorporate fluctuations of nucleon distri- 
butions. This shortcoming may result in favoring 
nuclei of smoothed simple shapes than in the real 
situation. 

(2) There is a strong possibility that some highly con- 
nected structures which have two or more substruc- 
tures in a period are neglected in these works be- 
cause only one structure is contained in a simula- 
tion box. 

It is not unnatural that the phases with highly con- 
nected nuclear and bubble regions are realized as the 
most energetically stable state [s^ [s^ . It is considered 
that, for example, a phase with perforated slablike nu- 
clei, which has negative x/V, could be more energetically 
stable than that with extremely thin slablike nuclei. The 
thin planar nucleus costs surface-surface energy which 
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FIG. 15: Euler characteristic density as a function of pth 
calculated for a; = 0.1 and p — 0.25po. The contribution of 
the nuclear surface can be observed as the plateau region. 

stems from the fact that nucleons bound in the nucleus 
feel its surfaces of both sides. The surface-surface energy 
brings about an extra energy increase in addition to the 
contribution of the surface energy. We have to examine 
the existence of the intermediate phases by more exten- 
sive simulations with larger nucleon numbers and with 
longer relaxation time scales in the future. 

2. Minkowski junctionals for X — Q.l 

In the case of a; = 0.1, the criterion for identification 
of the isodensity surface corresponding to the nuclear 
surface using the second derivative of V^(yOp, th)/^(Pp,th) 
does not work at higher densities. We thus use another 
method to calculate the normalized Minkowski function- 
als of the nuclear surface for x = 0.1. 

In Fig. El we have plotted the pth dependence of the 
Euler characteristic density x/V at p — 0.25po as an ex- 
ample. We can see that this curve consists of three com- 
ponents: the peaks of the lower pth region, the plateau 
region and the peaks of the higher pth region, which are 
due to dripped neutrons (thus these peaks cannot be ob- 
served for X — 0.5), nuclear surfaces and nucleons in nu- 
clei, respectively. These components can also be seen at 
the other values of p lower than However, we have to 
mention that the higher the density becomes, the smaller 
the plateau region gets, which means that the density 
contrast between the dripped neutron gas region and the 
nuclear matter region becomes obscure. Here, we take 
the mean values of the normalized Minkowski function- 
als in the plateau region as those for the nuclear surface, 
which are plotted as crosses in Fig. E| The error bars 
shown in this figure are the standard deviations of these 
quantities in the plateau region. Consistency between 
this method and the one using the second derivative of 
V{pp^th) /A{pp. th) has been confirmed for x = 0.3. 

Fig. El shows the resultant normalized Minkowski 
functionals for the nuclear surface at various values of 
p below Pin. The qualitative behaviors of u and (H) for 
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X = 0.1 are the same as those for x = 0.5 and 0.3; as p 
increases, u increases and {H) decreases (from positive to 
negative) almost monotonically in the density region of 
O-lpo ^ P ^ Pm- However, in the behaviors of A/V and 
x/V, quaUtative differences can be observed between the 
present case and the cases of x — 0.5 and 0.3. As can 
be seen in Fig. llGf b'). A/V increases almost linearly until 
just below Pm- 



The absence of the phases with cylindrical bubbles and 
with spherical bubbles in the phase diagram of x = 0.1 
(Fig. ISJ is well characterized by the behavior of x/V 
shown in Fig. Ilfif d'l. In the cases of x = 0.5 and 0.3, 
x/y increases from negative to positive with increasing 
density in the density region higher than that of the slab 
phase. However, for a; = 0.1, we cannot observe the ten- 
dency that x/V starts increasing even at just below p^. 
As a result, x/V keeps negative until matter becomes 
uniform. (Further discussion will be given in Section Ivl') 



Let us then consider the phase with slablike nuclei, 
which have not been obtained in the simulations for 
X = 0.1. If it were realized by using a longer relaxation 
time scale, it is expected to be obtained at p ~ 0.32 - 
0.34 po according to the behaviors of the {H) and x/V- 
However, we cannot see any signs from Fig. llSf b) that 
A/V stops increasing in this density region unlike the 
behaviors of A/V in the cases of x = 0.5 and 0.3. Here, 
we would like to mention that, according to the Landau- 
Peierls argument, thermal fluctuations are effective at de- 
stroying the long-range order of one-dimensional layered 
lattice of slablike nuclei rather than that of triangular 
lattice of rodlike nuclei and of bcc lattice of spherical 
ones. Thus, the melting temperature of the planar phase 
would be lower than the other phases, which leads to a 
longer time scale for formation of the slablike nuclei by 
the thermal diffusion. Therefore, a further investigation 
with a longer relaxation time scale is necessary to de- 
termine whether or not the phase with slablilkc nuclei 
is really prohibited in such neutron-rich matter in the 
present model. 



In Fig. El we have also plotted u and A/V for the 
density distribution p{r) as functions of pth like Figs. [T51 
and 1141 In comparison with Fig.^J the contribution of 
the dripped neutrons is shown more clearly in this case. 
We can see that the peak in the lower pth region due to 
the dripped neutrons combines to the peak in the higher 
Pth region. This behavior stems from the fact that a part 
of the dripped neutrons at lower densities are absorbed 
into nuclear matter region with increasing the density 
at fixed x; finally, all the neutrons are contained there 
in the uniform phase. We can also expect from the pth 
dependence of A/V that the phase with slablike nuclei 
might be obtained around 0.3poi where the slope of the 
plateau region is close to zero. 
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FIG. 16: Density dependence of the normalized Minkowski 
functionals for cold matter at a; = 0.1. 



PROPERTIES OF THE EFFECTIVE 
NUCLEAR INTERACTION 



Let us here examine the effective nuclear interaction 
used in this work. Structure of matter at subnuclear den- 
sities is affected by the properties of neutron-rich nuclei 
and of the pure neutron gas resulted from the nuclear 
interaction. Key quantities are the energy per nucleon 
e„ of the pure neutron matter, the proton chemical po- 
tential pL^^ in the pure neutron matter, and the nuclear 
surface tension i?surf- 

There is a tendency, especially in the case of neu- 
tron star matter, that the higher e„, the density pm at 
which matter becomes uniform is lowered. This is be- 
cause larger e„ tends to favor uniform nuclear matter 
without dripped neutron gas regions than mixed phases 
with dripped neutron gas regions. In the neutron star 
matter, there is also a tendency that the lower Pp"', the 

smaller p^. This is because —p^^ represents the degree 
to which the neutron gas outside the nuclei favors the 
presence of protons in itself. The quantity i^surf controls 
the size of the nuclei and bubbles, and hence the sum of 
the Coulomb and surface energies. With increasing i^surf 
and so this energy sum, pm gets lowered. 

It is important to check whether the effective nuclear 
force given by Eqs. l|B|l- l(T^ yields unrealistic values of 

these quantities or not. If e„, Ipp"'] and i?surf for the 
present model are unrealistically small in comparison 
with those for the other models, our results which have 
reproduced the "pasta" phases might be quite limited for 
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FIG. 17: Volume fraction u (upper panel) and surface area 
density A/V (lower panel) as functions of the threshold den- 
sity pth calculated for a: = 0.1 and various nucleon densities 
P- 



the present model Hamiltonian. 

In order to evaluate e„, we set 1372 neutrons in a sim- 
ulation box imposed of the periodic boundary condition. 
This system is cooled down by the QMD equations of 
motion with friction terms [see Eqs. H17(l ] until the tem- 
perature becomes ~ 1 keV. The resultant values of e„ arc 
plotted in Fig. ^1 We note that our results for p„ =0.2, 
0.6 and 1.0 po (the result for p„ ~ l.Opo is not plotted 
in Fig. I18|l coincide with the results for zero- proton ratio 
plotted in Fig. 9 of Ref. 

The values of e„ for the present model behave like 
those for the SkM Skyrme interaction especially in the 
density region of p„ < 0.13 fm~'^; they are close to the 
result of the variational chain summation obtained by 
Akmal, Pandharipande and Ravenhall [l^ at pn — po- 
The steep rise in e„ in the higher neutron density re- 
gion {pn <; 0.1 fm~^) compared to those obtained from 
the Hartree-Fock theory using various Skyrme interac- 
tions would help neutron-rich matter, which have larger 
dripped neutron density of pn ^ 0.1 fm^'^, to be uni- 
form. Therefore, we can say that this behavior of e„ 
for the present QMD model Hamiltonian suppresses the 
density region in which the "pasta" phases are the most 
energetically favorable in neutron star matter and in the 
case of x = 0.1 in the present study. We also note that 
e„ at lower neutron densities of p„ ^ 0.1 fm^'^ is rela- 
tively small. This, in turn, would lead to increase the 
density at which matter with lower dripped neutron den- 
sity (e.g., the case of a; = 0.3 in the present study) turns 
into uniform. 



. 18: The neutron density p„ dependence of the energy 
per nucleon e„ of the pure neutron matter. The solid squares 
show the result of the present QMD model Hamiltonian |2^. 
The dotted line denoted by SLy4 is the result from Ref. p4ll 
and the broken lines as marked by the other Skyrme inter- 
actions (FPS21, 1', FPS and SkM) are the results summa- 
rized by Pethick, Ravenhall and Lorenz The open stars 
and triangles denote the values obtained by Akmal, Pand- 
haripande and Ravenhall (3^ . and by Friedman and Pand- 
haripande [39^ . respectively. 



Next, we calculate the proton chemical potential p,)^' in 
the pure neutron matter. We use the cold neutron matter 
prepared for the above calculation of e„ as an initial con- 
dition. We insert a proton into this pure neutron matter, 
and then minimize the total energy by the frictional re- 
laxation method with fixing the positions and momenta 
of the other neutrons. The position of the inserted pro- 
ton is chosen randomly in the simulation box, and its 
momentum is chosen randomly from P < 30 MeV/c. We 

evaluate the difference in the total energy between 

that before the insertion of the proton and that after the 
optimization of the position and the momentum of the 
proton. 

In Fig. ^1 we plot p^'' for the present model Hamil- 
tonian. As can be seen from this figure, the result for 
the present model Hamiltonian generally reproduce the 
data of the other results obtained from the Hartree-Fock 
theory using the various Skyrme interactions at densities 
p ^ 0.1 fm"-^. At lower densities of p ^ 0.025 fm~'^, er- 
rors are quite large and data scatter significantly. This 
is because density fluctuations in pure neutron matter 
obtained by QMD would be unrealistically large at such 
low densities due to the fixed width of the wave packets 
in this model. However, it is noted that even in such a 
density region, our data are generally consistent with the 
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other results mentioned above. 
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FIG. 19: The neutron density dependence of the proton 
chemical potential fj'p^ in the pure neutron matter. The 
solid squares show the result of the present QMD Hamiltonian 
[2^ 1. The broken lines as marked by the Skyrme interactions 
(FPS21, r, FPS and SkM) are the results summarized by 
Pethick, Ravenhall and Lorenz |40||. and the solid line is the 
result of Sjoberg ^^l- The crosses denote the values obtained 
by Siemens and Pandharipande |4^. 

Finally, we turn to the surface tension, which af- 
fects energetically favorable nuclear shape most directly 
among the three quantities discussed here. We have cal- 
culated -Esurf from the energy change induced by the 
change in the area of planar nuclear matter. We use nu- 
cleon gas with proton fractions x — 0.1, 0.2, 0.3 and 0.5 
composed of 1372 nucleons. In the calculation of -Bsurf, 
the Coulomb interaction is excluded. 

In order to prepare slablike nuclear matter, we first 
cool down the above nucleon gas from U-qT ~ 20 MeV 
to ~ 0.2 MeV using Eqs. H17|) in a shallow trapping har- 
monic potential 



Viz) ^ 



(27) 



where k = 0.01 MeV fm~^. The simulation box here 
is imposed of the periodic boundary condition in the x- 
and y-directions, and is imposed of the open boundary 
condition in the z-direction. The box size L^.y in the x- 
and y-directions is set 20.26 fm. 

When the temperature reaches ^ 0.2 MeV, we remove 
the trapping potential and, except for the case with x = 
0.5, we change the boundary condition in the z-direction 
from open to periodic. The box size in the z-direction 
Lz is chosen so as to at least all nucleons dripped outside 
the nuclear matter region can be contained in the box: 



Lz = 92.32 fm {x = 0.1), 79.12 fm {x = 0.2), and 82.85 
fm [x — 0.3). After we relax the system for ~ 7000 
fm/c, we prepare three kinds of samples: one is nothing 
changed (sample 1) and the others have the area of the xy 
side of the simulation box increased (decreased) by 1 % 
with the total volume of the box kept constant [sample 2 
(sample 3)]. We further cool them down until fceT 0.1 
keV. The resultant nucleon density profiles for the sample 
1 of a; = 0.1 projected on the z-axis is shown in Fig. 
1201 As can be seen from this figure, Lz = 92.32 fm is 
much larger than the thickness of the slab d ~ 20 fm, 
and thus the volume of the dripped neutron gas region 
is almost the same among the three kinds of samples 
because the volume change in the nuclear matter region 
is negligible. It is also noted that d is much larger than 
the surface thickness dsurf — 5 fm, which ensures that 
the surfaces at both sides of the slab are separated well. 
Therefore, we can say that the energy difference between 
the sample 2 and 3 is just due to the difference in the 
surface area of the planar nuclear matter. Following the 
spirit of Ravenhall, Bennett and Pethick (hereafter RBP) 
pij l and by using the sample 1, we define the proton 
fraction Xi^ in the nuclear matter region as an averaged 
value for the region of the width of ~ 5 fm in the central 
part of the slab, where the proton and neutron density 
profiles ripple around constant values. 
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FIG. 20: The nucleon density profiles projected on the z-axis 
for the sample 1 of j: = 0.1. (The ^-coordinate is shifted.) 

We extract i?surf from the total energy E2 of the sample 
2 and of the sample 3 as follows: 



E2 — E3 
2{S2 ~ S3 



(28) 



where 5*2 (5*3) is the area of the xy side of the sample 2 
(sample 3) given by 20.26^ x 1.01 fm^ (20.26^ x 0.99 fm^). 
The factor 2 in the denominator represents, of course, 
the contribution of the two sides of the planar nuclear 
matter. As shown in Fig. [21 all the results plotted in this 
figure almost coincide with each other at Xi^ = 0.5, where 
uncertainty in the nuclear surface tension is rather small. 
These results deviate significantly at lower ccin, and the 
values of i?surf of the present calculation lie between those 
obtained by Baym, Bethe and Pethick (hereafter BBP) 



16 



and by RBP [H, which is based on the Hartree-Fock 
calculation using a Skyrme interaction, at Xi„ ~ 0.15 — 
0.35. Thus we can say that, in comparison with the result 
by RBP taken to be standard here, the contribution of 
Esuri of the present QMD model tends to favor uniform 
nuclear matter rather than the inhomogencous "pasta" 
phases for lower x-m of < 3.5. 
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FIG. 21: The nuclear surface energy per unit area (the surface 
tension) versus the proton fraction xin in the nuclear matter 
region. The solid squares are the values of the present QMD 
Hamiltonian 12^1 . the solid curve is the RBP result from their 
Hartree-Fock calculations |4^. and the dotted curve is the 
BBP result 0. 

Now let us discuss the fact that the intermediate phases 
with a negative value of x/^ s-^^g obtained instead of the 
bubble phases in a wide density region for x = 0.1. In 
considering this problem, we should note the general ten- 
dency that \fJ.^^\ increases and i?suif decreases as matter 
becomes neutron-rich (see Figs. I19land l21f) . The density 
dependence of A/V monotonically increases until pm as 
shown in Fig. 1161 would partly stem from the small value 
of i?surf ■ The apparently lower melting density pm in this 
case than in the cases of a; = 0.5 and 0.3, even though 
-Esurf is small, is due to a large which increases 

typically of order 10 MeV as p„ increases. The small 
-Esurf and the large l^p"^] in neutron-rich matter would 
help nuclear matter regions and neutron gas regions mix 
each other at the cost of small surface energy below pm. 
As a result, the structures with a negative x/V^ could 

be favored. According to Fig. E| the quantity fi^^ ob- 
tained for the present model Hamiltonian is consistent 
with those for the other Skyrme-Hartree-Fock calcula- 
tions in the relevant region of p„ ^ 0.08 fm~^. It is 
thus possible that a result which shows the structure of 
matter changes from negative x/^ to uniform without 



undergoing "swiss cheese" structure will be obtained by 
another calculation for neutron-rich matter using some 
framework without assuming nuclear shape. 

In closing this subsection, we summarize the conse- 
quences of the resultant e„, pp"^ and Esurf for the present 
QMD Hamiltonian. 

1. For symmetric matter (x — Xin — 0.5) 

According to Esurf at xi^ = 0.5, the present model 
is consistent with the other results, and is an ap- 
propriate effective interaction for the study of the 
"pasta" phases at x = 0.5. 

2. For neutron-rich cases 

For neutron-rich cases such as x ~ 0.1, in which the 
dripped neutron density p„ grows ^ 0.1 fm~^ just 
before matter turns into uniform, the present QMD 
model can be taken as a conservative one in repro- 
ducing the "pasta" phases. Its e„, fi^^ and E'surf 
act to suppress the density region of the "pasta" 
phases compared to other Skyrme-Hartree-Fock re- 
sults. 

3. For intermediate cases 

At intermediate proton fraction of x ~ 0.3, e„ of 
the present model acts to favor the inhomogencous 
"pasta" phases rather than the uniform phase and 
Esurf acts in the opposite way in comparison with 
other Skyrme-Hartree-Fock results. 

VI. ASTROPHYSICAL DISCUSSIONS 

Here we would like to discuss astrophysical conse- 
quences of our results. Pethick and Potekhin have 
pointed out that elastic properties of "pasta" phases with 
rodlike and slablike nuclei are similar to those of liquid 
crystals, which stems from the similarity in the geometri- 
cal structures 15j. It can also be said that the intermedi- 
ate phases observed in the present work are "spongelike" 
(or "rubberlike" for (H) < 0) phases because they have 
both highly connected nuclear and bubble regions shown 
as x/V < 0. The elastic properties of the spongelike in- 
termediate phases are qualitatively different from those of 
the liquid-crystal-like "pasta" phases because the former 
ones do not have any directions in which the restoring 
force does not act; while the latter ones have. Our results 
imply that the intermediate phases occupy a significant 
fraction of the density region in which nonspherical nu- 
clei can be seen (see Figs. |21 and EJ. According to Figs. 01 
and El we expect that the maximum elastic energy that 
can be stored in the neutron star crust and supernova 
inner core is higher than that in the case where all non- 
spherical nuclei have simple "pasta" structures. Besides, 
the cylinder and the slab phases, which are liquid-crystal- 
like, lie between the spongelike intermediate phases or the 
crystalline solidlike phase, and the releasing of the strain 
energy would, in consequence, concentrate in the domain 



17 



of these liquid-crystal-like phases. The above effects of 
the intermediate phases should be taken into account in 
considering the crust dynamics of starquakes and hydro- 
dynamics of the core collapse, etc. if these phases exist 
in neutron star matter and supernova matter. In the 
context of pulsar glitch phenomena, the effects of the 
spongelike nuclei on the pinning rate and the creep ve- 
locity of superfluid neutron vortices also have yet to be 
investigated. 

For neutrino cooling of neutron stars, some version of 
the direct URCA process which is suggested by Lorenz et 
al. Q, that this might be allowed in the "pasta" phases, 
would be suppressed in the intermediate phases. This is 
due to the fact that the proton spectrum at the Fermi 
surface is no longer continuous in the spongelike nuclei. 
An important topic which we would like to mention is 
about the effects of the intermediate phases on neutrino 
trapping in supernova cores. The nuclear parts connect 
over a wide region which is much larger than that char- 
acterized by the typical neutrino wave length ^ 20 fm. 
Thus the neutrino scattering processes are no longer co- 
herent in contrast to the case of the spherical nuclei, and 
this may, in consequence, reduce the diffusion time scale 
of neutrinos as in the case of "pasta" phases with simple 
structures. This reduction softens the supernova matter 
and would thus act to enhance the amount of the released 
gravitational energy. It would be interesting to estimate 
the neutrino opacity of the spongelike phases and the 
"pasta" phases. 

Finally, we would like to mention the thermal fluctua- 
tions with long wavelengths leading to displacements of 
"pasta" nuclei, which cannot be incorporated into the 
simulations using a finite-size box. Even if we succeed 
in reproducing the phase with slablike nuclei in neutron- 
rich matter in the future study, we should remind the 
above effect of thermal fluctuations to consider the real 
situation of matter in inner crusts of neutron stars. Fol- 
lowing the discussions in Refs. 0,0,113 by a hquid-drop 
model, it is likely that the extension of slablike nuclei is 
limited to a finite length scale of ~ 0(10^ — 10"^) fm in 
the temperature regions typical for neutron star crusts 
and supernova cores. 



Our results imply the existence of at least the phase with 
cylindrical nuclei in neutron star crusts because they cool 
down keeping the local thermal equilibrium after proto- 
neutron stars are formed and their cooling time scale, 
which is macroscopic one, is much larger than the relax- 
ation time scale of our simulations. 

In addition to these "pasta" phases with simple struc- 
tures, our results obtained here also suggest the existence 
of intermediate phases which have complicated nuclear 
shapes. We have systematically analyzed the structure 
of matter with two-point correlation functions and with 
morphological measures "Minkowski functionals" , and 
have demonstrated how structure changes with increas- 
ing density. Making use of a topological quantity called 
Euler characteristic, which is one of the Minkowski func- 
tionals, it has been found that the intermediate phases 
can be characterized as those with negative Euler char- 
acteristic. This means that the intermediate phases have 
"spongelike" (or "rubberlike" for (H) < 0) structures 
which have both highly connected nuclear and bubble 
regions. The elastic properties of the spongelikc interme- 
diate phases are qualitatively different from those of the 
liquid-crystal-like "pasta" phases. 

We have also investigated the properties of the effec- 
tive QMD interaction used in the present work in order 
to examine the validity of our results. Important quan- 
tities which affect the structure of matter are the energy 
per nucleon e„ of the pure neutron matter, the proton 
chemical potential fi^'^ in pure neutron matter and the 
nuclear surface tension i?surf- The resultant these quan- 
tities show that the present QMD interaction has gener- 
ally reasonable properties at subnuclear densities among 
other nuclear interactions. It is thus concluded that our 
results are not exceptional ones in terms of nuclear forces. 

Our results which suggest the existence of the highly 
connected intermediate phases as well as the simple 
"pasta" phases provide a vivid picture that matter in 
neutron star inner crusts and supernova inner cores has 
a variety of material phases. The stellar region which we 
have tried to understand throughout this paper is rela- 
tively tiny, but there has quite rich properties which stem 
from the fancy structures of dense matter. 



VII. SUMMARY AND CONCLUSION 

We have performed QMD simulations for matter with 
fixed proton fractions x = 0.5, 0.3 and 0.1 at various den- 
sities below the normal nuclear density. Our calculations 
without any assumptions on the nuclear shape demon- 
strate that the "pasta" phases with rodlike nuclei, with 
slablike nuclei, with cylindrical bubbles and with spher- 
ical bubbles can be formed dynamically from hot uni- 
form matter within the time scale of r ^ O(10'^ — lO'') 
fm/c in the proton-rich cases of a; = 0.5 and 0.3. We 
also demonstrate that the "pasta" phase with cylindrical 
nuclei can be formed dynamically within the time scale 
of r '--^ 0(10'^) fm/c for the neutron-rich case oi x = 0.1. 
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APPENDIX A: THE EWALD SUM FOR 
PARTICLES WITH A GAUSSIAN CHARGE 
DISTRIBUTION 

In calculating a long-range interaction, such as the 
Coulomb interaction, it is necessary to sum up all con- 
tributions of particles at a sufficiently far distance. The 
Ewald sum is a familiar technique for efficiently comput- 
ing the long-range contributions in a system with the pe- 
riodic boundary condition (see, e.g., Refs. [26ll45l| : recent 
mathematically careful discussion relating to the condi- 
tionally conv erge nce of the Coulomb energy can be seen, 
e.g., in Ref. |4g). The basic idea of the Ewald sum is 
that the contributions of particles in a long distance in 
real space can be calculated as contributions in the neigh- 
borhood in Fourier space: the contributions of particles 
in a short distance is summed up in real space and those 
of particles in a long distance is summed up in Fourier 
space. 



screening charge 
charged particle 



background charge 

compensating charge 




FIG. 22: Schematic picture of charge distribution in the 
Ewald sum for particles with a Gaussian charge distribution. 



Let us consider a system consisting of charged parti- 
cles which have a Gaussian charge distribution and a uni- 
form background charge which cancels the total charge of 
charged particles. These N particles are assumed to be 
in a cubic simulation box with volume V = i^ox '^hich 
is imposed by the periodic boundary condition. 

If every particle i with total charge Zi is surrounded by 
a Caussian charge distribution with total charge — Z^, the 



electrostatic interaction of particle i turns into a screened 
short-range interaction. Thus the total Coulomb energy 
^Coui of this system can be decomposed into as follows 
(see Fig. 1231: 



^Coul — ^short range + Ui — UscU , 



(Al) 



where Wshort range IS the sum of the Coulomb energy be- 
tween an each unscreened charged particle i and the other 
charged particles with screening charges, Ui is that be- 
tween an each unscreened charged particle i and com- 
pensating charges which cancels the screening Gaussian 
charges, and i4eif is the sum of spurious self interactions 
between an each charged particle i and its compensating 
charge i. 

Charge densities of the real charged particles p{r) and 
of the screening charges ps (r) can be written as 



n i—1 

N 



(A2) 



n i—1 



N 



n i—1 



N 



^^Ps,(n,*)(l-), 



(A3) 



n i=l 



where a and aEwaid are the reciprocals of the widths of 
charge distributions for a charged particle and a screening 
charge, respectively, and n denotes a position vector of a 
periodic image normalized by Lbox- Thus a distribution 
of screened charges Pscreened is 



/Osereened(r) = p(r) -I" Ps{r) 



(A4) 



The electrostatic potential 0short range (i") due to 

^screened 

(r) can be obtained as 



J 



N , 

Pshort range 

n V— 1 ^ 



erfc(ya |r - (r,; -I- Lboxn)|) erfc(^Q;Ewaid |r - (rj + Lboxn)|) 



|r - (r^ iboxn)| 



|r - (r^ Lboxn) 



(A5) 



r 



because a solution of the Poisson equation 



IS 



1 (P- 

--r^{r(l){r)) = 47r 
r ar'^ 



(r) = Z^I^(^+ const., 



(A6) 



(A7) 



where erf(x) and erfc(x) are the error function and the 
complementary error function, respectively, and they are 
defined as erf(a:) = /p^exp(— s^) ds and erfc(a;) = 

1 — erf(x). Here, we determine the constant C so that 
the average value of 0short range in the simulation box V 
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be zero: 



Thus, 0short range (r) leads to 



X 



^short range (^) ^ ^ 



a aEwald 



(A8) 
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Pshort range 



n i=l 
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crfc(^rt |r - (r,: + L],^^n)\) crfc(^nK«-a.id |i' - (r, + Ii,oxn)|) 
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a Q^Ewald 



I Pavr ) 



where the average charge density Pavr of the charged par- 
ticles is defined as 



'y 

i=l 



(AlO) 



(A9) 



r 



The total Coulomb energy Wshort range between an each 
unscreened real charged particle i and the other charged 
particles with a screening charge can be calculated as 



short range 



lf]|/rf\/dVp(„=o,i)(r) 

i=l I 



n' J 



erfc(i/a |r - (r^ + Lboxn)|) erfc(y'Q;Ewaid |r - (fj + -^boxn)|) ' 



+ Z/boxn) 



+ 



|r - {vj + Lboxn) 



+ 



a Q!Ewald / 
9 EE Ir.- - 



ZiZj 



|ri - (r^- + iboxn)| 



erfc ( ^ / " '^^'^'^'d I J,. _ (j,^. _|_ Lboxn)| 



a + aEwald 



-erfc (rj +Lboxn)| 



V / TT 



TT \ 



O- Q^Ewald , 



r 



(All) 



where the primes on the summations mean the terms Fourier transforming the charge distribution pi (r) yields 
i = j at n = are excluded. 



Next, we calculate Ui , which is the sum of the Coulomb 

energy between an unscreened charged particle i and a 
charge density pi(r) which consists of the compensating 
charges and the background charge: 



= y E ^ie-*--^e-'^i;^ - pavr (A13) 



Pl(r) = - Ps(r) - Pavr 
AT 



Using pi(k) and the Poisson equation (V^0i(r) 
-47rpi(r)) in the Fourier form 



n j=l 
Pavr ' 



fcVi(k)=47rpi(k) , (A14) 
(A12) we can at once obtain the electrostatic potential (pi due 
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to the charge density pi. 
Mr) = 5]0i(k)e^'^"- 



N 

k^O 3=1 

The term with k = is canceled due to charge neutrahty. 
Thus the Coulomb energy lAi is given by 



1 ^ r 



a\3/2 



-a|r-ri|' 



Mr) 



2 ^ ^ yfc2 

k#0 i,] 



I ik-(ri-r,) 



X e 



(A16) 



We have to subtract a sum of spurious self interactions 
^sGif,i between a charged particle i and its compensating 
charge from Ui. According to Eq. (|A7I) . an electrostatic 
potential i^Causs due to a Gaussian compensating charge 
is 0Gauss(r) = Zi erf(^aEwaid r)/r, thus the self interac- 
tion of particle i reads to 



Wself,i = y" Z, (^^) ^ e '"'VGauss(r-) d?T 
Vtt/ V a + aEwald 



and hence, 



Wself — - 



N 



2 ^ 

1=1 



self ,i 



a QfEwald 
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ttEwald 



- N 

E 



zf 



(A17) 



(A18) 



Finally, the total Coulomb energy can be calculated by 
Eqs. (EH, jIlT|l . (|Il6|) and llllgl) . The positive back- 
ground charge does not appear explicitly because the av- 
erage value of (pshort range within the simulation box is set 
to be zero. 

The total Coulomb energy per particle ecoui and the x 
component of the Coulomb force fi^x acting on a particle 
for various values of aEwaid are plotted in Fig. In 
this calculation, we use 1024 positive charged particles 
(protons) distributed randomly in a simulation box of 
ibox = 39.59 fm (i.e. pp = O.lpo), which is imposed 
of the periodic boundary condition. Figs. [231 (a) and [231 
(b) show the results for point charges and for Gaussian 
charge distributions, respectively, which are calculated 
for the same configuration of the particle positions {r^}. 
The width of the Gaussian charge distributions is set to 
be a = 1/2L with L = 2.1 fm^, which corresponds to the 



width of the wave packets in the QMD model used in this 
work. 

We note that there are plateau regions of ecoui and fi ^ 
whose values do not depend on aEwaid- These constant 
values give results to be obtained. We note that the range 



-0.2 



-0.4 



g 0.05 

o 

v_ 



-0.05 
-0.1 

-0.15 
-0.2 




(a) point particles 

_i 1 1 I 1 < < < L 



"T 



(b) Gaussian =7 

(a=l/2L, L = 2.1) 1<™=10 

k^.x=14 

1 1 I 1 < < < \ 1 1 1 1 I 1 1 i_ 



0.5 1 1.5 



FIG. 23: The total Coulomb energy per particle ecoui (in 
units of MeV) and the x component of the Coulomb force 
(in units of MeV/fm) acting on a particle obtained from 
Eqs. lAlL IIAllL IIAI6I 1 and IAI8II as a function of .^OEwaW. 
We use 1024 protons distributed randomly in a box of Lbox = 
39.59 fm (i.e. pp = O.lpo). The results shown in (a) and (b) 
are calculated for the same configuration of particle positions 
{vi} but for different values of the width a of the Gaussian 
distribution; (a) a = (point charge) and (b) a — 1/2L with 
L = 2.1 fm^ 



of -/ctEwaid of the plateau regions become larger with in- 
creasing fcmax, where fcmax is the cutoff radius in the unit 
of 27r/Lbox for the summation in Fourier space. As can 
be seen from Fig. [231 ctEwaid dependences of ecoui and 
fi^x for the present QMD model with a finite width of 
the Gaussian charge distributions are weaker than those 
for the point charges. These features are also confirmed 
for different proton number densities of 0.2 and 0.3 pq. In 
our simulations, aEwaid is set 13 or 14, which are consid- 
ered to be large enough to calculate the total Coulomb 
energy per particle in accuracy less than 0(1) keV, which 
is the typical value of the energy difference between suc- 
cessive "pasta" phases in neutron star matter obtained 
by previous works. 
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